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ABSTRACT 


A  new  method  of  numerical  solution  for  time-harmonic 
field  scattering  problems  involving  inhomogeneous  bodies  of 
revolution  is  presented.  A  simple  single-loop  multiport 
feedback  system  is  driven  by  the  incident  field.  The  forward 
path  transfer  matrix,  found  by  use  of  the  finite-element 
method,  translates  the  total  near-field  into  the  body-surface 
field.  The  feedback  path  matrix  is  established  by  the  magne¬ 
tic  vector  potential  formulation*,  it  translates  the  body- 
surface  currents  into  the  scattered  near-field.  The  thin- 
wire  scattering  problem  is  formulated  in  terms  of  this  field 
feedback  method.  Results  of  this  computation  are  compared  in 
graphical  form  to  those  found  by  Hallen's  integral  equation 
solution  of  the  scattering  problem  for  several  cases  of 
interest . 
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I.  INTRODUCTION 


The  objective  of  this  endeavor  is  to  investigate  the  per¬ 
formance  of  a  new,  and  potentially  quite  powerful,  technique 
for  the  computation  of  electromagnetic  scattering  problems. 

For  reasons  that  will  later  become  apparent,  this  new  technique 

3 

will  be  denoted  the  "field  feedback  formulation"  (FFF  or  F  ) . 

•  3 

The  work  reported  here  is  a  first  application  of  F  ,  where 
the  calculation  of  scattering  from  straight  thin-wires  has 
been  performed.  This  is  an  almost  classical  electromagnetic 
field  problem  that  has  been  routinely  solved  by  use  of  inte¬ 
gral  equations  for  many  years.  Such  a  cononical  problem  was 

3 

chosen  as  a  testing  vehicle  for  the  newly  conceived  F 
analysis  method  since  alternate  integral  equation  solutions 
are  readily  available  for  comparisons. 

The  formulation  developed  here  employs  a  combination  of 
the  finite  element  numerical  method  and  the  classical  magnetic 
vector  potential.  The  finite  element  method  (FEM)  is  an 
application  of  variational  techniques  which  implements  a 
weighted  residual  approximation.  The  high  accuracy  of  the 
finite  element  method  is  used  to  advantage  in  reducing  the 
amount  of  computing  required,  as  compared  to  other  numerical 
methods.  Use  is  also  made  of  coupled  azimuthal  potentials 
(CAP).  The  CAP  formulation  restricts  the  class  of  approachable 


problems  to  those  which  are  axially  symmetric,  but  does  allow 
the  treatment  of  material  inhomogeneities  and  multiwavelength 
dimensions . 

Thin-wire  assumptions  are  made  which  greatly  simplify 
the  current-scattering  problem.  These  assumptions,  however, 
are  by  no  means  critical  to  the  solution  method.  The  formu¬ 
lation  should  be  readily  extendable  to  a  large  class  of  impor¬ 
tant  interior,  radiation,  and  scattering  problems. 


II.  THE  PROBLEM  AND  METHOD  OF  SOLUTION 


Electromagnetic  scattering  problems  usually  involve 
a  known  incident  field,  such  as  a  uniform  plane  wave.  The 
scattered  field  can  be  found  from  the  current  induced  by  the 
presence  of  the  incident  field,  but  the  induced  current  can 
be  found  only  after  the  total  field  is  known.  Since  the 
total  field  is  just  the  sum  of  the  incident  and  scattered 
fields,  a  definite  difficulty  has  been  encountered.  Know¬ 
ledge  of  the  induced  current  is  needed  to  find  the  scattered 
field,  while  at  the  same  time,  the  scattered  field  must  be 
known  to  determine  the  induced  current . 

This  is  clearly  seen  in  the  form  of  Maxwell's  equations 
for  the  scattering  problem.  Assume  that  a  scatterer  with 
properties  Y^=c^+ jwe^ ,  Z^=jwy,  is  in  free  space  where 
Yo  =  jwuo,  ZQ  =  j(jJuo,  and  in  the  presence  of  sources  J Q ,  Mq. 

Then  the  curl  equations  have  the  form 

VxH  =  7  +Y,F  (la) 

o  1 

VxE  =  -M  -A.  H  (lb) 

o  1 

where  the  induced  electric  and  magnetic  currents  are 

T,  =  (Y. -Y  )E  (2a) 

l  1  o 

FL  =  (Z.-z  )H 
1  1  o 
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■hm 


and  the  scattered  magnetic  field  is 


H 


scat 


=  FT  -  H 


incident 


(3) 


The  electric  field  relationship  is,  of  course,  similar. 

The  simple  and  compact  nature  of  the  simultaneous  equa¬ 
tions  (1)  and  (2)  is  highly  appealing,  but  of  little  practi¬ 
cal  use  in  obtaining  numerical  answers  to  a  problem  because 
of  the  excessive  computing  required.  Yet,  while  the  require¬ 
ment  to  have  two  unknowns  cannot  be  escaped,  the  equations 
relating  these  unknowns  can  at  least  be  rendered  more 
tractable. 


A.  THE  FIELD  FEEDBACK  FORMULATION 

As  will  be  shown  later,  two  simple  matrix  equations  can 
be  derived  which  relate  the  magnetic  and  electric  fields  at 
a  number  of  points  along  the  surface  of  the  scatterer  to  the 
fields  at  points  along  a  boundary  some  distance  away.  These 
equations  are  independent  and  work  in  opposite  directions. 
That  is,  from  the  first  equation,  given  the  total  boundary 
field  values,  e_t  and  h  ,  the  scatterer  surface  total  field 
values,  e  and  h  ,  can  be  found: 

— OJ  —il) 

[e Jh  ]T  =  V  Ce.  I  h.  3T  (4a) 

— — 0)  ~  — t  j  — t 

T 

The  column  vectors,  [e.|h.]  ,  contain  the  field  values  at  the 

i 

f  T 

computation  points.  The  second  equation  takes  Ce^lh^]  and 

T 

yields  the  scattered  field  values  Ces|hs]  along  the  boundary 


(4b) 


CeJh  ]T 

s  i  s 


T  [e.  ! h, ,] 


— w  -w' 


Equations  (4a),  (4b),  and  (3)  comprise  a  simple  feedback 
system,  as  in  Figure  1. 


Figure  1.  Field  Feedback 

What  remains,  then,  is  to  construct  the  matrices  V  and  T 
to  fit  the  problem  at  hand.  It  should  be  pointed  out  that 
the  magnetic  field  alone  is  involved  in  the  wire-scatterer 
problem  as  developed  here.  So  [e  !h  ]  is  replaced  by  h  . 

The  field  translation  matrices  V  and  X>  though  developed  in 
jxactly  the  same  way,  will  be  considerably  simplified  as 
compared  to  the  general  scatterer  case. 

B.  THE  WIRE  SCATTERING  PROBLEM 

Given  a  perfectly  conducting  thin  wire  in  free  space,  and 
a  plane  wave  of  known  amplitude  incident  from  a  known  direction, 
it  is  required  to  find  the  current  induced  on  the  wire  and  the 
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field  scattered  by  the  wire  at  any  point  in  space.  Finding 
the  induced  current  is  equivalent  to  finding  the  field  on  the 

^  A  _  A 

wire  surface,  since  J  =  n  x  H  ,  where  n  is  the  unit  vector 
normal  to  the  surface.  Normalized  coordinates  will  be  used 
throughout  the  problem,  and  are  given  by 

(R,Z,<}>)  =  (K  r,  Kqz,  <p). 

The  coordinates  (r,Z,i)  are  standard  circular  cylindri¬ 
cal  coordinates,  and  K  =  2 tt / X  is  the  free  space  wavenumber  of 
the  time-harmonic  field.  The  situation  is  depicted  in 
Figure  2.  Use  of  these  normalized  coordinates  will  allow 
the  formation  of  a  compact  and  useful  form  of  the  CAP  equations. 
The  wire  radius  is  a,  given  in  terms  of  XQ,  the  free  space 

wavelength.  Then  the  normalized  radius  is  R  =  2iTa/X  .  The 

a  ~  o 

wire  length  is  L,  also  given  in  wavelengths.  So  Zt=2ttL/X  . 

The  incidence  angle  of  the  plane  wave  is  a,  measured  as 
shown  in  Figure  2.  Since  perfect  conductivity  is  assumed, 
fields  and  currents  will  exist  only  on  the  surface  of  the 
wire.  And  since  a  very  thin  wire  will  be  assumed,  where  a 
<X  /10,  the  current  at  the  wire  ends  may  be  set  to  zero. 

The  problem  will  not  be  solved  in  terms  of  the  usual 
electric  and  magnetic  fields,  E  and  H,  or  the  scalar  poten¬ 
tial,  V.  Instead,  use  will  be  made  of  the  magnetic  vector 
potential,  A.  The  current  and  scattered  fields  will  be 
presented  and  calculated  in  terms  of  the  CAP  formulation 
discussed  below. 


C.  COUPLED  AZIMUTHAL  POTENTIALS 

The  CAP  formulation  [3]  has  been  developed  to  handle 
previously  intractable  time-harmonic  field  problems  involving 
continuously  and  discretely  inhomogeneous  bodies  of  revolu¬ 
tion.  With  this  limitation  of  axial  symmetry,  all  electro¬ 
magnetic  field  components  can  be  expanded  via  exponential 
Fourier  series  in  the  <p  coordinate.  The  model  fields  can 
then  be  represented  by  two  coupled  azimuthal  potentials. 

These  potentials  satisfy  a  system  of  coupled  partial  differ¬ 
ential  equations,  as  well  as  a  variational  criterion,  which 
will  be  developed  here. 

Using  the  normalized  coordinates  in  Figure  2,  the  total 
fields  are  decomposed  into  azimuthal  modes: 
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T 


m 


E(R,Z , 4? ) 


(R,Z)exp(  jm4>) 


noH(R,Z,<J>) 


h  (R,Z)exp(  j  nx4> ) 


(5a) 


(5b) 


where  no  =  120ir.  By  substitution  of  these  expansions  into 
Maxwell’s  equations,  it  can  be  shown  that  all  modal  field 
components  can  be  generated  from  two  coupled  azimuthal 
potentials  i^1(R,Z,m)  and  !//2(3,Z,m),  using 

*N  /V 

<(>  x  em  =  jfni(m<t)x7ii)^-Ru  VilJj )  (6a) 


♦  '  em  =  ^i/R  (6b) 

A  A 

<p  x  hm  =  jfn(m<pxV!l/2->-RerV'ip1)  (6c) 

A 

<p  *  hm  =  ^/R  ( 6d) 


A 

where  <j>  is  the  unit  azimuthal  vector,  7  is  the  two  dimen¬ 
sional  gradient  operator,  and 


fm[R,Z]=[ur(R,Z)er(R,Z)R2-m2]“1  (7) 

The  Euler-Lagrange  variational  technique  will  be  dis¬ 
cussed  more  fully  in  the  exposition  of  the  finite  element 
method  below.  It  is  sufficient  to  note  here  that  this  tech¬ 
nique  finds  the  solution  of  the  system  partial  differential 
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equations  as  the  stationary  point  of  a  functional  of  the 
potentials  and  their  first  derivatives.  This  stationary 
point  is  defined  as  the  point  (^o’  ^20^  which  the 
variation  of  F,  5F,  vanishes  for  perturbations  (5i/>^,  5  \f/^) 
about  that  point. 

The  functional  F  is  a  surface  integral  of  the  Lagrangian 
L,  where 

A  /v 

L  =  (R£r7^1+Tn<})x7ijj2)  +  V^2  •  ] 

-  (v^2+v^22)/R  (8) 

The  functional  F  is  integrated  over  the  constant  azimuth 
planar  meridian  cross  section  Q  of  the  volume  in  which  the  solu¬ 
tion  is  of  interest. 

It  is  convenient  at  this  point  to  change  variables,  and 
introduce  a  simplifying  assumption.  From  (6b)  and  (6d)  we 
have  e^m(R ,Z) s^^/R  and  h^m(R,Z)  =  iJ^/R.  So  finding  only  the 
azimuthal  components  of  the  modal  fields,  e^m  and  h^,  will 
yield  a  complete  solution  through  (6a)  and  (6c).  Now,  since 
a  thin  wire  is  assumed,  the  phase  variation  of  the  incident 
field  from  side  to  side  on  the  wire  may  be  neglected.  Then 

A 

the  induced  current  will  be  a  function  of  Z  alone;  J^^Z. 
Therefore  the  scattered  fields  will  be  azimuthally  invariant, 
so  that  E(R,Z,<J>)  =  eQ(R,Z)  and  nQH(R,Z  ,4> )  =hQ(R,Z  ) .  Further- 

A 

more,  since  J^nxf^H^Z  and  since 
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1  *  HJ, 


*0' 


4)0 


The  number  of  modes 


jue  jwe 

it  follows  that  e,  =0,  and  h  =  h 

4>o  ’ 

has  now  been  reduced  to  one,  the  h,  mode. 

4>o 

Observe  that  (6a)  and  (6c)  no  longer  involve  coupled 
potentials.  This  allows  the  Lagrangian  to  be  decoupled,  and 
it  may  now  be  written 


L  =  f  {V(Rh.n)-CRir,7(Rh.n)]}-y  Rhf  . 
o  <po  pr  <po  r  <po 


(9) 


So  the  functional  is 


=  /  L<R.z.V>7V)dR<iz 

n 


(10) 


It  remains  to  find  h^o(R,Z)  to  make  F  stationary  in  ft  for 
given  Dirichlet  boundary  conditions .  This  will  be  done  by 
use  of  the  finite  element  method. 
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III.  COMPUTATIONAL  PROCEDURES 


The  constructions  of  the  V  and  T  matrices  have  little 
in  common  other  than  the  coordinate  system.  The  finite 
element  method  is  used  to  generate  V.  T  is  derived  from 
the  magnetic  vector  potential  in  this  perfect  conductor  case, 
while  if  a  dielectric  body  were  used,  both  electric  and 
magnetic  vector  potentials  would  be  required  to  construct  T. 

The  use  of  numerical  methods  obviously  implies  that  the 
current  and  fields  will  be  found  only  at  discrete  points. 

As  discussed  in  Section  II. C.,  only  points  in  a  constant 
azimuth  meridian  plane  need  be  considered.  The  region  between 
the  wire  surface  and  the  boundary  at  which  the  fields  are  to 
be  calculated  is  divided  into  triangular  elements.  These  are 
areas  of  integrarion  in  the  finite  element  method.  The  divi¬ 
sions  and  the  nodal  numbering  scheme  are  shown  in  Figure  3. 

There  are  a  total  of  points  along  the  wire,  and  Ng  points 
along  the  boundary. 

A.  FINITE  ELEMENT  METHOD 

This  method  is  one  whereby  the  region  of  interest  is  divided 
into  a  number  of  overlapping  subregions  of  a  specified  shape. 
Then  the  problem  solution  function  9  is  approximated  by 
defining  the  solution  as  a  linear  combination  of  the  products 
of  unknown  coefficients  with  "basis"  functions  defined  over 


over  each  of  the  overlapping  subregions.  These  basis 
functions  are  usually  sectionally  linear,  and,  in  very  simple 
terms,  are  used  to  "fill  in"  the  approximation  between  points. 
A  weighted  residual  approximation  of  the  solution  9,  commonly 
known  as  the  Ritz  method,  can  then  be  made  by  the  application 
of  the  discretized  Euler-Lagrange  variational  procedure, 
subject  to  the  appropriate  boundary  conditions.  The  rigorous 
development  of  each  of  these  parts  of  the  finite  element 
method  is  not  necessary  here,  but  the  background  required  for 
a  basic  understanding  of  the  method  will  be  given  below. 


Figure  3.  Solution  Points 
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Basis  Function  Expansions 


In  order  to  solve  a  problem  numerically,  some  method 
of  approximating  the  solution  using  a  finite  number  of 
unknowns  must  be  found.  Also,  the  method's  convenience  must 
be  balanced  against  the  desired  accuracy.  A  general  term 
for  such  an  approximation  is  a  "basis  function  expansion" 
of  a  function  0  [2].  The  type  of  expansion  used  in  the 
finite  element  method  for  a  one-dimensicnal  problem  is  given 
by 

N 

9NCX)  =  S  9iu.(X)  (11) 

i  =  l 

where  9^  is  the  discrete  version  of  9,  the  9^  are  the 
values  of  9(X)  at  the  discrete  points  X^ ,  and  the  IMX) 
are  the  basis  functions  at  the  points  X^.  Linear  basis 
functions  are  usually  employed  in  FEX  calculations,  though 
quadratic  or  higher-order  polynomial  functions  could  be 
used  to  obtain  greater  accuracy.  Linear  basis  functions  for 
a  one-dimensional  problem  are  shown  in  Figure  4.  Note  that 
the  node  spacing  is  irregular.  This  allows  tailoring  of  the 
expansion  to  better  approximate  the  solution,  and,  in  two 
or  higher  dimensional  cases,  may  be  necessary  because  of  the 
shape  of  the  solution  region.  If  neither  of  these  are 
necessary,  regular  node  spacing  can  be  used  with  some  advan¬ 
tage  to  reduce  complexity.  The  basis  fucntions  u^  are  seen 


1 


to  have  unit  value  at  their  respective  nodes  x^,  and  are 
zero  for  x^_^>x>x^+^  >  lj^i^M  • 


Y  f 

1 


X 


Figure  4.  Linear  Basis  Functions 


A  two-dimensional  linear  basis  function  often  used  is 
the  "pyramid"  function;  its  two-dimensional  projection  is 
shown  in  Figure  5.  A  rectangular  region  2  in  the  x-y  plane 
has  been  discretized  by  division  into  the  triangular  elements 
shown.  The  boundary  3232  in  2  of  the  function  u3  2(x,y)  has 
been  outlined  to  make  the  shape  clear. 

The  function  has  a  value  of  one  at  the  center  node,  and 
is  zero  along  and  outside  the  periphery.  Its  equation  can 
be  written 


u32(x,y) 


<x>y)  , 


j 


a,f 


(12) 


where  the  component  functions  are  planar  segments  of  the  form 


Figure  5.  Two-Dimensional  Finite  Element  Mesh 

Then  a,3,Y  must  be  functions  of  (x,y)  so  that  u^ 2 (x3  ,y2 )  =  1 , 

u32  ^xi  ,yj  ^  =  0  ’  3*2,  fr°r  integer  (i,j). 

As  in  the  one-dimensional  case,  the  basis  functions  over¬ 
lap;  that  is,  a  function  is  centered  around  each  node  in  the 
region  or  on  the  boundary.  Of  course,  on  the  boundary  of  71 
at  least  half  of  the  function  is  "lost"  outside  71.  The 
overall  appearance  of  the  basis  function  array  in  three  di¬ 
mensions  might  be  imagined  as  a  crystalline  surface  of  multi¬ 
faceted  peaks  centered  over  the  nodes,  with  raised  valleys 
running  between.  If  a  function  9(x,y)  were  defined  over  71, 
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its  expansion  could  be  written 


7  6 

v*’y)  *  H  2 

i=j  5=1 


(14) 


If  the  exact  values  of  0-.  were  known,  9V  would  be  a  multi- 

i  j  N 

faceted  approximation  to  8^  ^  ,  equal  only  at  the  nodes  x^ . . 
Then,  within  the  limitations  of  machine  accuracy  and  computing 
time,  the  greater  the  node  density  the  better  the  approxima¬ 
tion  8^  will  be. 

The  wire  scattering  problem  expands  the  field,  h^Q,  by 

use  of  these  pyramid  basis  functions.  Referring  to  Figure  3, 
it  is  seen  that  the  region  H  covered  by  the  triangular  ele¬ 
ments  will  not  contain  any  complete  pyramids,  but  this  is  not 
really  necessary.  It  will  prove  convenient  to  use  a  sequen¬ 
tial  node  coordinate  system  rather  than  the  (i,j)  Cartesian 
description  of  (12).  Let  n  be  the  node  coordinate,  where 
l_<n_<N ,  M=N^+Ng.  Then  the  approximating  expansion  h^  of  h^Q 
can  be  written 


N 

h  ,.(R,Z)  =  f  h  u  ( R ,  Z ) 
<pN  J  n  n 

n=l 


(15) 


Note  that  the  top  and  bottom  wire  nodes  are  not  represented. 
This  is  because  of  the  assumption  that  h^o=0  at  these  points. 

Two  field  vectors  can  now  be  defined.  The  nodal  wire- 
surface  field  is 


h  =  h 
— w  — n 


n  =  1,  N 


w 


(16) 
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The  nodal  boundary  field  is 


-B  "  -n  ’  n  V1*  N 


These  vectors  contain  the  approximate  values  of  h,  at  the 

vpO 

nodes  of  the  triangular  mesh. 

2  .  Euler-Lagrange  Variational  Technique 

Consider  the  differential  equation  A9  =  f,  which  together 
with  certain  boundary  conditions,  describes  some  physical  pro¬ 
blem  in  a  region  Q.  The  complex  scalar  wave  (Helmholz) 

2  2 

equation,  A  =[7  +k  ]9=Q  is  an  example.  The  differential  equa¬ 
tion  is  related  to  the  functional  F(9),  where 


F(  9  ) 


=  /  L(9,  D^Q ,  D  9,  x,  y)dH 


This  relationship  is  a  fundamental  result  of  the  calculus  of 
variations,  and  simply  stated,  is  the  fact  that  the  functional 
F(9)  is  stationary  at  the  "point"  9q  which  also  satisfies  the 
differential  equation.  That  is,  if  the  particular  function  9q 
can  be  found  which  causes  the  derivative  (or  first  variation) 
of  F(9)  to  vanish,  then  F(9)  is  said  to  be  stationary,  and 
9q  is  said  to  be  the  stationary  point  of  F(9);  9q  is  then 
the  solution  of  the  equation  A9=f  [5]. 

Observe  that  the  vanishing  of  dF/d9  implies  that  a 
minimum,  maximum,  or  saddle  point  of  F  occurs  9  =  9q.  in  many 
physical  problems  F  is  associated  with  the  potential  energy 
of  systems  and  thus  finding  9q  is  sometimes  termed  the 
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*  *■ 


"minimization"  of  F.  The  problem  now  presents  a  choice  of 
computational  methods.  The  differential  equation  might  be 
approximated  by  a  discrete  system,  say  finite  differences, 
or  the  variational  integral,  F,  can  be  "minimized"  over  a 
finite  number  of  trial  functions,  as  is  done  by  the  finite 
element  method. 

It  is  necessary  in  this  procedure  to  find  the  Lagrangian 
L  for  use  in  F ( 9 )  [4].  This  L  must  be  such  that  its  insertion 
in  the  Euler-Lagrange  equation  produces  the  differential 
equation  that  is  to  be  solved.  As  an  example,  the  Euler- 
Lagrange  equation  for  the  two  independent  variable  case  is 
given  by 


3L  3  r  3L  i  3  r  3L  ■, 

30  “  3X  L  ao  eJ  ■  3Y  L  SD  SJ 
x  y 


(19) 


There  is  no  definite  way  to  discover  L  for  any  given  problem, 

but  it  usually  yields  to  informed  guessing. 

The  problem's  boundary  conditions  present  no  unmanageable 

difficulties  in  this  formulation.  If  0  is  specified  on  the 

boundary  (Dirichlet  boundary  conditions)  nothing  else  need 

be  done.  If  9  is  unspecified  on  the  boundard,  the  stationary 

solution  will  automatically  fulfill  the  Neumann  boundary 

condition,  =0,  also  known  as  the  natural  boundary  condition, 
dn 

Boundary  conditions  of  a  more  complicated  type  will  require 
the  functional  F  to  be  modified,  but  this  procedure  will  not 


23 


be  necessary  here.  It  can  be  easily  shown  by  using  the  field 
generating  equations  in  [6],  that  the  natural  boundary  condi¬ 
tion  here  results  in  a  zero  tangential  electric  field  on  the 
wire  surface. 

3 .  Application  to  the  Thin-Wire  Problem 

As  discussed  in  Section  II,  the  driving  force  in  the 

scattering  problem  is  the  incident  plane  wave.  Then,  the 

boundary  conditions  to  be  applied  are  the  natural  condition 

along  the  wire,  for  nodes  n=l,  N  ,  and  a  Dirichlet  condition 

along  the  outer  boundary  of  Q;  h.Q=ht  (total  field).  The 

Dirichlet  condition  h,  =0  has  also  been  imposed  on  the  too 

$o  r 

and  bottom  wire  nodes.  These  correspond  to  boundary  condi¬ 
tions  along  the  top  and  bottom  node  rows  of  a  larger  mesh, 
as  in  Figure  5 . 

The  minimization  of  F  for  the  wire  scattering  problem 

is  carried  out  as  follows.  F  is  discretized  by  substituting 

h,.r  for  h.  .  which  yields  the  aDDroximate  functional 
<pN  bo  J 

F(h*N)  =  /  [fofV(RS,'()'R“rV<RV>)-R'Jrh*N]dRdZ  l20) 

The  exact  functional  F(h^Q)  is  made  stationary  about  the 
solution  h^Q  by  differentiating  F  with  respect  to  h^Q.  This 
makes  F  stationary  about  every  point  on  h^.  Since  the  appro¬ 
ximate  solution  is  to  be  found  only  at  a  finite  number  (N^?) 
of  points,  the  functional  needs  to  be  made  stationary  only  at 
those  points.  This  is  done  by  a  simultaneous  solution  of  the 
Nw  stationary  condition  equations 
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25 


where 


[J?Ru.*VRu  -Ru.u  JdRdZ  (24b) 

R  1  n  i  n 

The  matrix  equation  to  be  used  for  computation  can  be 

formed  directly  from  (24).  First,  recall  that  the  basis 

functions  u  are  zero  on  and  outside  of  the  boundary  . 

n  y  n 

This  means  that  the  product  of  non-overlapping  basis  func¬ 
tions  is  zero.  For  example,  in  Figure  3,  u~*u, =0.  Thus 

^  _L 

u . u  =0  ,  I n-i I  >  2  (25a) 

and 

VRu. -VRu  =0  ,  I n-i 1  >2  (25b) 

i  n  1  1  — 


This  behavior  is  very  useful  in  that  it  produces  banded 
coefficient  matrices.  Special  techniques  can  be  used  to 
invert  large  banded  matrices  very  quickly.  Writing  out  a 
few  equations  should  make  the  matrix  construction  clear. 


i  =  1 : 

Cllhl+C12h2  + 

i=  2  : 

C21hl+C22h2+C 

i  =  3 : 

C32h2+C33h3+C 

C1N  +1  hN  +1  +  C1N  +2hN  +2  =  0 

w  w  w  w 

23h3  +  C2N  +2hN  +2+C2N  +3hN  +3  = 
w  w  w  w 

34h4  +  C3N  +3hN  +3+C3N  +4hN  +4  = 
w  w  w  w 


0 

0 

(26) 


r 


r 


The  system  equation  can  therefore  be  written 

Ah  =  -3hD  (27) 

w  B 


where  the  N  xN  coefficient  matrix  A  contains  the  basis 
w  w 

function  integrals  for  i  =  l,N  ,n=l,M  and  the  N  xN  +1  matrix 
°  w  5  w  w  w 

B  contains  the  integrals  for  i=l,N  ,n=N  +1,N.  Both  A  and  B 
^  w  w  - 

are  diagonal  and  have  the  forms 


X  X 

X  X 

B: 

XXX  0 

XX  0 

XXX 

X  X 

0  ’x  X  X 

0  XX 

X  X 

X  X 

hw  due  to  the  boundary  field  h3  is  now 

(28) 

Equation  (23)  is  the  forward  path  in  the  system  of  Figure  1. 

The  actual  calculation  of  the  A.  ,  B.  will  be  discussed  in 

m  in 

Section  III.C. 

B.  SCATTERED  FIELD  INTEGRAL 

1 .  Vector  Potential  Formulation 

The  development  of  the  feedback  path  is  less  compli¬ 
cated  than  the  forward  path  construction  given  in  III. A.  The 


The  wire  field 
known,  since 

h  =-A-1BhD  =  V  h_ 
— w  ~  B  -  — t 
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thin-wire  approximation  is  applied  to  the  familiar  classical 
vector  potential  K,  and  a  matrix  equation  formulation  built 
on  a  one-dimensional  basis  function  expansion.  The  result 
will  be  an  approximation  giving  the  near  field  produced  by  a 
wire-surface  field  or  current. 

From  the  inhomogeneous  Helmholz  equation 

V2A  +  k2A  =  -7  (29) 


it  is  seen,  since  J=JZZ  by  assumption,  that  A-A^Z .  Then  it 
follows  from  the  vector  potential  defining  equation 

H  =  V  x  A  (30) 


that 

a*  vr-z)  <3i) 

The  azimuthal  invariance  of  H  has  been  made  explicit  in  (31). 

A  commonly  used  expression  for  A  [1]  is  given  as  an 
integration  of  a  weighted  point-source  Green's  function  G(r' ,r) 
over  the  current  source  where 


G(r ' ,r ) 


_1_  e-jko!r“rJ 

4n  [F-r 1 | 


Then 


A(r) 


J(r ' )G(r'  ,r)dV ' 


This  general  situation  is  shown  in  Figure  6. 


(32) 


(33) 
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Figure  6 


Scattered  Field  Matrix  Equation 


The  formulation  is  easily  applied  to  the  thin-wire 
case.  Recalling  that  h^  =  nQH  ,  and  assuming  that  the 
current  J„  lies  along  the  Z  axis,  A,  for  the  thin  wire,  is 


n0Az<R,Z) 


/./v 


e-JlDl 

4tt|D| 


R  d<pdZ ' 


2  2 

where,  |  D  |  =  +(Z-Z)  ,  as  seen  in  Figure  7,  and  the  integra 
tion  is  over  S,  the  wire  surface.  Now,  putting  (34)  in  (31) 
yields 


1/2ZT 


(R,Z)=- 


(Z’  )[- 


|D| 


This  h^Q(R,Z)  is  approximately  the  scattered 

field,  h  (R,Z),  radiated  by  the  wire  carrying  an  induced 

current  J„  .  This  wire  current,  in  the  form  of  h,  C Z ’ )  ,  must 
L  (pO 

be  approximated  by  a  basis  function  expansion  so  that  hg(R,Z) 
can  be  computed.  Since  h^CZ’)  is  a  function  of  Z  only,  the 
one-dimensional  linear  basis  functions  shown  in  Figure  4  may 
be  used.  Then,  for  the  same  nodes  along  the  wire,  as  in 
Section  III .A. , 


(Z»  ) 


h  u  (Z '  ) 
n  n 


(36) 


Notice  that  (36)  defines  the  same  nodal  surface-field  vector 
as  (16).  Substitution  of  (36)  into  (35)  produces 
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R  R 

hg(R,Z)=-|- 


»„  f 

t  h"/ 

n=i  J  Z' 


n+l 


u  (Z’)[  — U-  +  — ^]e~^  '^dZ' 

|D|  IDT 


n 


n-1 


(37) 


The  limits  of  integration  reduce  from  +ZT  to  Z*  , 

—  1j/  n— ■  1 5 


Z*  , ,  because  the  basis  function  u  (Z1)  is  zero  for  Z’  ,>Z'>Z  ,, 
n+l  n  n-1—  —  n+l 


Equation  (37)  can  be  used  to  construct  the  T  matrix  by  noting 
that  (37)  gives  the  scattered  field  at  one  point  (R,Z).  Then, 
by  choosing  a  number  of  points  in  some  region,  a  discrete 
field  approximation  can  be  constructed.  The  points  chosen 
will  be,  unsurprisingly,  just  those  Ng  points  along  the  finite 
element  mesh  border.  So  the  matrix  equation  is 


*s 


=  T  h 


-  — w 


(38a) 


where 


R  R 
a 


'bn 


/  n+l 
Z '  n-1 


u  (Z  '  )  [ — + 

n  ID? 


r]  e"j  lDUz' 


(38b) 


is  the  coefficient  of  the  n-th  value  of  the  wire  field,  and 

T.  h  is  one  term  of  the  b-th  value  of  the  scattered  field, 
bn  n 

Finally,  since  h^  =  hg  +  h^nc,  an  expression  for  the 
total  field  in  terms  of  the  incident  and  wire  fields  is 


Ht  =  — inc 


♦  T  h 
-  — w 


(39) 
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The  scattered  field  h  found  in  (38)  is  a  near  field 

— s 

expression.  It  is  valid  at  any  distance  from  the  wire,  but 
this  behavior  does  not  extend  to  the  entire  problem.  If  hg 
were  to  be  found  at  some  large  distance  from  the  wire,  simply 
extending  the  radial  width  of  would  not  suffice.  This  is 
because  FEM  calculation  accuracy  degrades  as  the  distance 
from  the  wire  to  the  boundary,  BT  ,  increases.  The  magnitudes 
of  the  electrostatic  and  induction  fields  decrease  rapidly 
as  R  increases,  so  the  point  is  quickly  reached  where  they 
become  effectively  zero  as  far  as  computer  calculations  are 
concerned.  The  resulting  incorrect  total  field  values  are  a 
part  of  the  FEM  R-related  inaccuracies.  The  greater  part  of 
the  FEM  error  is  found  in  the  basis  function  expansion  of  h^ . 

The  planar-section  basis  functions  used  are  unable  to 
accurately  model  the  near  field  over  a  large  radial  distance, 
where  "large",  in  this  case,  may  be  1/100  of  the  wavelength. 

To  solve  the  far-field  problem,  the  near-field  solu¬ 
tion  must  first  be  obtained  for  the  wire  field  h  .  Then  (38) 

~"W 

may  be  used  to  find  hg  along  any  boundary.  For  typical  far- 
field  problems,  assumptions  may  be  made  to  reduce  the  integral’s 
complexity.  Let  8  be  the  angle  of  the  far-field  point, 
measured  from  the  Z  axis  at  the  origin.  Then  R=kQr  sin8  , 
where  r  is  the  spherical  radius  of  the  point.  The  far-field 
assumption  is 

| D I  =  kQr  -  Z’  cos3  (in  phase  term  of  (38b))  (40a) 

| D |  =  kQr  (in  denominators) 


(40b) 


So  the  radiation  field  along  a  spherical  arc  of  radius  r  is 
given  by  (38),  where  T^n  is 


u  (Z’)e^Z’COs3 


dZ  ' 


(41) 


C.  PROBLEM  IMPLEMENTATION 

Combining  equations  (27)  and  (39)  gives 

h  =  - ( A+ST ) _1B  h.  „  (42) 

This  is,  of  course,  the  same  equation  as  is  represented  in 
Figure  1,  where  V  =  -A  1B.  The  system  "transfer  matrix" 
is  defined  by 

h  =  (I  -  V  •  T)"1  •  V  •  h.  (43a) 

— w  ~  —me 

Simple  matrix  algebra  leads  to  a  simplification  of  this 
result 

h  =  -(A  +  B  •  T)”1  *  B  •  h.  (43b) 

— w  ~  ~  —me 

Multiplying  this  equation  by  (A+BT)  yields  (42),  so  the  system 

equation  as  derived  does  behave  as  a  feedback  system. 

The  only  variable  remaining  to  be  identified  is  h.  ,  the 

—me 

incident  field  vector.  The  incident  plane  wave  h.  can  be 

1  me 
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expressed  in  the  normalized  cylindrical  coordinate  system  as 


h^(R,Z)  =  je"^Zcosa  JQ ' (Rsina)  (44) 

where  J  '  is  the  derivative  of  the  zeroth-order  Bessel  func- 
o 

tion.  Since  R  will  be  kept  small,  as  discussed  in  Section 
III.B.,  (44)  may  be  simplified  by  replacing  .J  1  (Rsina) 
with  the  first  few  terms  of  its  series  expansion 


J  ' (Rsina) 
o 


,Rsina  ,x ,  Rsina,  3  ,  Rsina, 
v  2  K  16  384  ’ 


(45) 


Then  h.  may  be  calculated  at  the  boundary  nodes  by  use  of 
(44)  and  (45). 

Now  that  expressions  for  ail  of  the  quantities  in  (42)  have 

been  found,  the  computer  program  can  be  devised  to  carry  out 

the  calcuations .  The  use  of  an  existing  routine,  CSMINV ,  for 

the  inversion  of  (A+3T)  reduces  the  number  of  major  computations 

to  two:  equation  (24b)  for  the  A  and  B  elements,  and  equation 

(38b)  for  the  scattered  field  coefficients  T,  . 

bn 

The  integration  of  (34b)  becomes  algebraicly  messy,  but  is 
otherwide  straightforward.  Because  of  the  basis  function’s 
behavior,  the  integration  area  is  reduced  from  !)  to  322 ,  the 
region  of  non-zero  overlap  of  u^  and  u  .  This  integration 
over  322  may  then  be  written  as  a  sum  of  integrations  over  the 
triangular  elements  comprising  322,  since  is  piecewise  planar 
in  the  elements.  Also,  note  that  C.  is  independent  of  Z  in 


the  sense  that  the  integration  for  any  pair  (i,n)  of  nodes, 
will  he  the  same  if  that  pair  is  translated  along  the  wire. 

This  means  that  one  wire  node  may  be  selected,  and  the  inte¬ 
grals  concerning  that  node  and  its  immediate  neighbors  will 
be  the  only  ones  required  to  construct  A  and  B.  For  example, 
select  i=2.  Then  C  anc^  <"2M  +3  are  t^ie  on^-^ 

integrals  that  need  to  be  calculated.  (Mode  three  is  not 
being  neglected.  Equation  (24b)  shows  that  C_,  =  C.~,  but 

this  is  ^ust  C,,  moved  up  one  element.  So  C-,=C„..) 

2  1  r  2  3  2  u 

In  this  simple  narrow  mesh  case,  the  integration  of  (24b) 
is  required  over  at  most  three  elements  (C .,,,),  and  the  size 
and  shape  of  all  elements  is  the  sane.  For  many  problems 
this  will  not  be  the  case.  For  large  meshes,  and  especially 
for  skewed  elements,  it  becomes  very  convenient  computa¬ 
tionally  to  change  from  a  global  to  a  local  node  system. 

Instead  of  moving  from  note  to  node,  calculating  all  the 
coefficients  about  each  node  one  at  a  time,  an  algorithm 
is  devised  to  move  from  element  to  element,  finding  the 
integrals  involving  the  three  corner  nodes  over  each  element. 
This  procedure  effects  a  considerable  simplification  in  the 
"bookkeeping"  necessary  to  keep  track  of  the  nodal  indices 
and  coordinates. 

The  local  node  system  employed  here  is  shown  in  Figure  8. 
The  integrals  of  the  combinations  of  local  nodes  (k,l)  are 
found  over  a  upper  and  lower  element.  Then  the  integrals 
can  be  placed  in  their  proper  global  positions.  As  an  example, 
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Upper 

Element 


Lower 

Element 


Figure  8.  Local  Node  Indices 

consider  the  upper  element  with  global  nodes  (2,3,N  +3).  The 

integrations  of  (24b)  over  U  are  found  for  local  node  pairs 

(1,1) ,(1,2)  ,(1,3) ,(3,1)  ,(3,2)  ,  and  (3,3).  Then  the  local 

node  pair  (1,1)  integral,  /(1,1)  is  put  in  A  as  one  part  in 

three  of  C22,  /( 1,2)  is  one  part  in  two  of  C2N  +  3  in  B ,  and 

w 

/(1,3)=C23  in  A.  Also,  /(3,1)=C32,  /( 3,2)  is  part  of  C3,j  +3, 

W 

and  /(3,3)  is  part  of  C33.  The  remaining  three  local  node 
integrals  over  U  are  not  needed  in  this  case,  but  would  be 
used  in  larger  meshes. 

For  the  thin  mesh  used  here,  the  local  node  procedure 
improves  only  slightly  the  algorithm  efficiency  as  compared 
to  a  global  node  procedure.  The  principle  reason  for  its  use 
was  that  an  existing  and  successful  subroutine  was  available 
to  do  local  node  integrations.  This  routine  is  named  VAP^LA 
in  the  program.  It  requires  only  the  local  node  ordering 
scheme  and  the  coordinate  (R,Z)  of  each  local  node  to  produce 
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a  three  by  three  matrix  containing  the  integrations.  As 
discussed  above,  the  wire-scattering  problem  needs  only  two 
calls  to  VARELA  to  obtain  the  necessary  values.  The  sub¬ 
routine  ABLOAD  then  fills  A  and  B  with  these  integrals. 

The  calculation  of  the  coefficients  T,  is  done  by  the 

bm  J 

familiar  Gaussian  quadrature  method,  using  eight  unequally 
spaced  points.  The  interval  of  integration  in  (38b)  is  split 
into  two  segments,  as  in  Figure  9.  Then  a  change  of  variable 
is  affected  so  that  the  limits  of  integration  in  each  sub¬ 
interval  are  X=+l  to  -1.  The  Gaussian  procedure  then  gives 


+  1 

L  f 


(x)dx  ^  1  w . [ f ( x . ) +  f ( -x . ) ] 

/  -i  it  l 

i  =  1 


where  the  values  +v.  and  the  weight  factors  w.  are  tabulated 
in  several  references  such  as  [6].  The  code  used  to  imple¬ 
ment  these  integrations  is  contained  in  subroutine  GFINT . 

Equation  (42)  is  coded  in  two  subroutines,  CFORM  and 
ZCUPR.  Routine  CFORM  produces  the  matrix  product 
(A+BT)-  3=C ,  while  ZCURR  calculates  the  incident  field 


vector  h.  according  to  (44)  and  (45),  then  finds  the  wire 
—me  J 

field  vector  h  =C  h. 

~w  ~  -  me 


U  (Z) 
v  n 


Jn-1  n  n+1 

Figure  9 


IV.  RESULTS  AND  ANALYSIS 


In  many  other  time-harmonic  field  problems,  FEM  calcula¬ 
tions  yield  good  results  with  perhaps  2Q  to  30  nodes  per 
wavelength.  With  higher  node  densities,  say  50-75  perwaveler.gth, 
the  solution  breaks  down  due  to  roundoff  error  in  the  compu¬ 
tations.  As  will  be  shown,  this  does  not  appear  to  be  the 
case  with  the  field  feedback  formulation.  Nodal  densities 
as  high  as  100  per  wavelength  have  been  used  without  diffi¬ 
culty.  Indeed,  these  high  densities  are  required  in  order  to 
achieve  good  convergence  of  the  formulation.  As  yet,  no 
satisfactory  answer  for  this  behavior  has  been  found. 

The  only  other  parameter  available  to  affect  the  conver¬ 
gence  of  the  solution  is  ST  ,  the  wire-to-boundary  distance. 

As  mentioned  in  section  III. 3. 2.,  this  parameter  can  be  ex¬ 
pected  to  be  critical  in  achieving  a  good  field  approximation. 
Computed  results  show  this  to  be  true. 

A  great  deal  of  effort  was  expended  in  looking  for  errors 
in  the  formulation  before  the  solution  behavior  with  respect 
to  and  N^  was  fully  appreciated.  It  is  still  possible,  no 
matter  how  unlikely,  that  some  error  still  exists.  If  an 
error  is  present  in  the  formulation,  or  code,  it’s  effect  is 
quite  peculiar.  By  comparison  with  other  FEM  problems,  the 
possible  error  causes  two  to  three  times  the  "normal"  node 
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density  to  be  needed,  but  otherwise  seems  to  have  no  effect. 
That  is,  the  results  are  satisfactory  in  all  respects  save 
node  density  required  for  convergence. 

The  primary  comparison  case  chosen  was  a  one-half  wave¬ 
length  wire  with  a  radius  of  1/50  wavelength.  This  case  was 
a  compromise  forced  by  excessive  required  computing  for  long-^ 
wires  and  questionable  integral  equation  comparison  results 
with  shorter  wires.  These  other  cases  were  looked  at,  and 
some  example  results  are  included  here.  Their  investigation, 
however,  has  not  been  exhaustively  completed  due  to  lack  of 
time . 

A.  INTEGRAL  EQUATION  COMPARISON 

The  basis  used  for  judging  the  convergence  of  the  FFM 
computation  was  a  standard  and  successful  type  of  integral 
equation  formulation,  Hallen's  equation  for  scattering  by 
a  thin  conductive  cylinder.  This  method  assumes  a  Z  directed 
current  density,  and  finds  this  wire  current  by  a  simultaneou 
solution  of  two  expressions  for  the  vector  potential,  A.  One 
equation  is  just  equation  (33).  The  second  is  a  different 
form  of  equation  (20) 

V2A  +  k2A  =  j we  Es(r)  (47) 

_ s  _  ... 

where  E  (r)  is  the  scattered  electric  field.  Simple  manipu¬ 
lations  yield  the  integral  equation  for  I(z'),  [ 7 J 
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+  CsinkZ+D  cos  kZ  -  U(Z)  (48) 

where  A  is  an  impedance  per  unit  length,  C  and  D  are  ccn- 
starts  of  integration,  and  U  contains  the  incident  field. 
Equation  (43)  is  coded  in  subroutine  INTEQ . 

Results  for  the  1/2  x  1/50  wavelength  case  are  shown  on 

3 

pages  A-l  through  A-38  of  Appendix  A.  F  calculations  are 
shown  as  a  solid  line,  while  integral  equation  computations 
are  displayed  as  the  "X"  curve.  The  parameters  of  each  case 
are  given  below  the  graph.  The  magnitude  plots  are  scaled 
in  milliamperes  along  x,  and  wavelengths  along  y.  The  wire 
is  to  be  seen  lying  along  the  y  axis,  and  is  conveniently 
demarcated  by  the  zero  end  values  of  the  integral  equation 
curve.  The  current  phase  plots  have  the  same  y  axis  units, 
and  are  scaled  in  degrees  along  x. 

Graphs  A-l  through  A-24  show  convergence  with  an 
increasing  number  of  nodes.  Values  of  are  20,40,70  and 
100,  with  three  plana  wave  incidence  angles  (ALPHA)  used  in 
each  case.  In  these  plots,  the  plane  wave  is  incident  from 
the  upper  right  for  ALPHA=30  and  60  degrees.  ALPHA=90°  is 
broadside  incidence.  Convergence  of  the  current  magnitude 
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is  plainly  seen  at  all  incidence  angles.  For  a=90  ,  the 

error  is  15.5%  at  the  wire  center  for  N  =20,  but  only  2.5%  for 

w 

N  =70.  Similarly,  the  chase  error  decreases  from  13°  tc  1°. 
w 

Comparing  the  N^=  7 0  and  Nw=100  cases,  it  is  seen  that  the 
extra  nodes  make  very  little  difference  in  the  quality  of  the 
solution.  Convergence  has  been  essentially  accomplished  at 
Nw=70,  and  the  additional  memory  and  computing  time  for  Nw=100 
are  not  really  necessary.  The  convergence  is  less  satisfying 
for  incidence  angles  other  than  90°.  At  ALPHA=30°,  for  example, 
the  magnitude  is  not  skewed  quite  enough,  while  the  phase  error 
increases  toward  the  wire  bottom.  This  situation  does  not 
imorove  as  N  increases. 

W 

Graphs  A-25  through  A-38  show  FFM  convergence  as  a  function 

3 

of  B.  .  All  F  calculations  are  made  with  N  =28.  The  behavior 

Li  W 

shown  in  these  plots  is  characteristic  of  all  cases  of  wire 
length  and  radius.  For  much  too  large,  the  magnitude  is 
too  low  and  the  phase  is  off  by  as  much  as  90°.  As  3,  decreases, 
the  magnitude  rises,  overshoots,  and  then  settles  back  towards 
the  correct,  convergent  value.  The  phase  moves  smoothly  down 

3 

to  the  correct  value.  Notice  that  the  F  curves  are  nearly 
identical  for  DR=DZ/32  and  DR=DZ/64.  This  implies  that 
convergence  has  been  reached  by  the  time  DR  has  decreased  to 
DZ/32,  or  B.=.0005X  .  For  this  case  the  convergence  point  of 
B.  is  actually  about  BT=.001X  ,  as  shown  in  A-l  through  A-25. 
Results  for  several  other  cases  of  wire  length  and  radius 
are  shown  in  Appendix  B.  The  one  and  two  wavelength  wire 
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calculations  show  that  more  nodes  are  required,  assuming  that 

B^=.001Xq  is  still  adequate  since  has  not  changed.  If  70 

nodes  are  required  for  a  good  answer  at  Z^=Xq/2,  then  for 

two  wavelengths  as  many  as  300  nodes  might  be  necessary. 

Similarly,  well  over  100  nodes  could  be  needed  for  the  one 

wavelength  wire.  The  one  wavelength  case  indicates  that 

3 

this  is  probably  true.  The  F  solution  for  Nw=100  is  percep¬ 
tibly  better  than  the  N  =70  solution,  esoecially  for  off-axis 
incidence  angles . 

Appendix  C  contains  graphs  for  two  different  cases: 

A^  =  l/100Xo  while  =  A q / 2  or  Z^=Aq/4.  Graphs  C-l  through  C-4 

again  show  convergence  value  is  no  longer  B^.0011^,  but  is 

smaller  than  3  =.00051  . 

L  o 

Plots  C-7  through  C-10  show  the  quarter-wavelength  wire 
case.  It  appears  from  the  FFM  curves  that  convergence  has 
been  almost  reached  for  N^  =  4  0 .  This  calls  into  question 
the  accuracy  of  the  integral  equation  solution.  This  possible 
inaccuracy  should  not  be  too  surprising,  however,  since  the 
wire  is  relatively  thick  for  it's  length,  and  the  integral 
equation  does  much  better  for  very  thin  wires. 

3 

An  "error  isolation"  check  was  made  on  the  F  system  by 
testing  the  feedback  path.  A  fictitious  sinusoidal  wire 
current  was  generated  and  its  radiated  field  calculated  via 
h  =T  h  .  This  field  was  compared  to  that  given  in  reference 
[33.  The  comparison  is  shown  in  plots  C-ll  and  C-12.  The 
two  field  calculations  are  equal,  to  the  fourth  significant 
digit .  This  indicates  that  the  feedback  path  formulation  and 
program  coding  is  correct. 
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B.  ERROR  ANALYSIS 


Consideration  of  the  near  field  behavior  reveals  the  reason 
for  the  critical  dependence  of  the  solution  on  3^.  The  magne¬ 
tic  field  magnitude  in  the  very  near  field  is  dominated  by  a 
2 

1/R  term,  as  can  be  deduced  from  equation  (37).  This  is 

depicted  in  Figure  10.  Also  shown  are  two  dashed  lines  which 

are  the  linear  basis  function  expansion  of  the  field.  It  is 

clear  from  the  figure  that  the  approximation  of  the  field 

deteriorates  as  A^  decreases,  if  is  held  constant.  Thus 

it  is  seen  why  3T=.001Ao  is  adequate  for  A^=.C2Aq  while 

3. =.00051  to  obtain  the  best  convergence  with  A=.01A  .  It 
L  o  L  o 

is  noted,  however,  that  even  at  this  best  convergence  value 

3 

of  3^,  the  F  solution  is  not  as  close  to  the  integral  equation 

solution  as  in  the  A=.02A  ,  B=.001A  case.  This  fact  can 

L  o  L  o 

also  be  explained  by  consideration  of  the  basis  function 

approximation's  behavior.  While  decreasing  B^  improves  the 

field  approximation  very  near  the  wire,  the  approximation 

becomes  worse  in  the  neighborhood  of  R=B^* 

It  thus  appears  that  there  is  a  limit  to  the  usefulness 
.  3 

of  this  formulation  of  the  F  solution  to  the  wire-scattering 

problem.  The  primary  difficulty  lies  in  the  inaccuracy  of 

the  linear  basis  function  expansion  used  in  the  FEM  calculation. 

A  much  better  basis  function  may  contain  a  term  proportional 
2 

to  1/R  .  This  would  greatly  reduce  the  sensitivity  of  the 
L’ 


\ 


solution  to  B 


There  is  another  source  of  error  which  is  inherent  in 
any  sort  of  computer  solution:  roundoff.  Simple  calculations 
will  be  affected  in  the  sixth  or  seventh  significant  digit  by 
the  use  of  single-precision  arithmetic  on  the  IBM-360.  If 
only  a  few  operations  were  required,  this  error  could  be 
ignored.  But  when  solving  large  systems  of  equations  many 
thousands  of  operations  are  required,  and  accuracy  can  be 
severely  degraded. 

A  rough  guide  to  the  degree  of  trust  to  be  placed  in  a 
solution  can  be  had  in  the  condition  number  of  an  inverted 
matrix.  This  number  provides  an  upper  bound  on  error  propa¬ 
gation  in  a  system  due  to  the  inaccuracy  of  finite-precision 
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arithmetic.  For  example,  in  the  system  X  =  A  ‘'"Y,  a  condition 
number  of  1000  for  A~^  would  cause  errors  in  the  third  signi¬ 
ficant  digit  of  Y  to  appear  as  errors  in  the  third,  second  or 
even  first  significant  digit  of  X. 

This  condition  number  is  a  measure  of  the  similarity  of 
the  system  equations.  In  the  wire-scatterer  case,  as  the 
solution  points  are  placed  closer  together,  the  field  values 
at  these  nodes  become  more  indistinguishable.  Hence  the 
equations  relating  the  field  values  at  neighboring  nodes 
become  more  nearly  the  same.  The  system  solution,  then, 
becomes  sensitive  to  errors  in  the  less  significant  digits. 

Condition  numbers  of  system  matrices  for  the  cases  dis¬ 
cussed  here  range  from  about  250  to  1000.  Cases  with  fewer 
nodes  and  larger  values  of  3r  have  lower  condition  numbers , 
while  as  N  increases  and/or  8.  decreases,  the  condition 

W  L 

number  increases  rapidly.  This  is  another  reason  why  a  good 
solution  is  more  difficult  to  achieve  for  very  thin  or  long 
wires .  The  use  of  a  basis  function  more  suited  to  the  field 
behavior  would  allow  larger  values  of  to  be  employed.  This 
would  lower  the  system  condition  number  significantly,  and 
improve  the  solution's  accuracy  and  trustworthiness.  Other 
than  this  and  keeping  Nw  as  small  as  possible,  there  is 
nothing  that  can  be  done  to  avoid  the  problems  reflected  in 
the  condition  number. 

Even  if  the  errors  mentioned  above  could  be  eliminated, 

3 

the  results  of  both  the  F  and  integral  equation  calculations 
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would  still  not  represent  the  fields  with  complete  accuracy. 

The  assumptions  made  in  both  formulations  are  not  strictly 
true  for  any  wire  of  finite  thickness.  Current  flows  on 
the  surface  of  the  wire,  not  in  the  center,  as  assumed  in 
the  vector  potential  formulation. 

Field  phase  variation  from  side  to  side  of  the  wire  is  not 
zero,  so  there  is  an  azimuthal  component  of  current  as  well  as  the 

A  A 

Z-directed  current  assumed.  And  because  this  ^-directed 
current  exists,  the  current  is  not  zero  at  the  wire  ends. 

But,  since  essentially  the  same  types  of  approximations  are 
made  for  both  solution  methods,  it  is  not  possible  to  ascer¬ 
tain,  from  the  comparisons  made  here,  the  degree  to  which 
either  solution  represents  the  actual  fields. 
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V.  CONCLUSIONS 


As  a  result  of  this  investigation,  the  field-feedback 

formulation  has  been  shown  to  be  a  sound  technique  for  the 

computation  of  electromagnetic  scattering  problems.  It  has 

3 

been  demonstrated  that  F  calculations  compare  very  favorably 
with  results  obtained  by  the  conventional  integral  equation 
formulation.  The  success  of  this  initial  development  is  not 
unqualified,  however.  The  two  deficiencies  previously  dis¬ 
cussed,  basis  function  unsuitability  and  excessive  number  of 
required  nodes,  could  present  formidable  difficulties  in 
problems  where  a  priori  knowledge  of  field  behavior  is  not 
available,  or  where  bodies  large  compared  to  a  wavelength  are 
used . 

The  discrepancy  between  required  and  expected  node 

density  is  particularly  troublesome,  as  this  fact  is  as  yet 

unexplained.  If  high  node  densities  are  actually  required 

by  the  method  and  are  not  just  the  result  of  some  formulational 

error,  then  the  attendant  large  matrices  will  limit  the  use- 
3 

fulness  of  F  computations  because  of  the  great  amount  of  com¬ 
puter  memory  required. 

In  any  case,  the  flexibility  of  the  finite-element  method 
when  used  in  the  coupled-azimuthal  potential  formulation 
should  allow  numerous  heretofore  difficult  problems  to  be 
attacked.  The  ease  with  which  mixed  boundary  conditions  can 


be  met  for  complicated  inhomogeneous  scattering  bodies,  and 


the  accuracy  and  simplicity  of  the  vector  potential  formu¬ 
lation  promise  a  wide  applicability  for  the  field  feedback 
formulation. 
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ZL=0. 5: AL=0. 02: BL=0. 001 : ALF fl  =  30 : NW=  I  00 

A- 2  3 


71 


73 
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-020  -015 


-OLD  -aos 


OQO  005 


4- 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

"oTo  +ojs' 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


L 


X-SCRLE=5.00E*01  UNITS  INCH. 

Y-SCRLE=5. 00E-01  UNITS  INCH. 

CURRENT  PHASE  :FEM/T=SOLID  :INT.EO=PLUS 
ZL  = . 5  :  HL= 1 / 1 00  :  DR=DZ/4  :  RLFA=90 

A-  2  6 


74 


-025  -020  -015  -010  -005  000  000  010  019  020 


1 


X-SCfilE=5. 00E-04  UNITS  INCH. 

T-SCRLE=5.  CJQE-0  1  UNITS  INCH. 

CURRENT  MAGNITUDE  :FEM/T=SGLID  :INT.EQ=X 
ZL  = . 5  :  AL= 1 / 1 00  :  DR=DZ/9.8  :  RLF A  =  90 

A-  2  7 
75 


I 


on  -020  -013  -010  -003 


X-SCRLE  =  5.00E-Oti  UNITS  INCH. 
Y-SCALE=5.00E-01  UNITS  INCH. 


CURRENT  MAGNITUDE  :FEM/T  =  SOLID  :  I  NT . EQ  =  X 
ZL  = .  5  :  AL=l/50  :  DRA=[?gZ/2  :  ALFA  =  90 


77 


X-SCRLE=5. 00E+01  UNITS  INCH. 

T-SCflLE  =  5.  OOE-0 1  UNITS  INCH. 

CURRENT  PHASE  :FEM/T=SGLID  :INT.EQ=PLUS 
ZL=. 5  :  AL=l/50  :  DR=DZ/2  :  ALFA=90 

A-  3  0 
78 


X-SCRIE=5 . 00E+0 1  UNITS  INCH. 
Y-SCRLE=5.00E-0!  UNITS  INCH. 


CURRENT  PHASE  :FEM/T=S0LID 


ZL  = .  5  :  AL=l/50  :  DR  =  DZ/4 

A-  3  2 


I  NT . EQ  =  PLUS 
RLFA  =  90 


80 


X-SCRLE=5. 00E+0 1  UNITS  INCH. 

T-SCAIE=5, OOE-O 1  UNITS  INCH. 

CURRENT  PHRSE  :FEM/T=50LID  :INT.EG=PLUS 
ZL  = .  5  :  RL=l/50  :  DR=DZ/8  :  RLFR=90 

A- 3  4 


32 


-ota 


010 


-005 


+ 

+ 

+ 

+ 

+ 


X-SCRLE=5. 00E+0 1  UNITS  INCH. 
Y-SCRIE=5. OOE-O 1  UNITS  INCH. 


CURRENT  PHASE  :FEM/T=S0LID 
ZL  = .  5  :  RL=l/50  :  DR=DZ/64 

A-  3  8 


I  NT . EQ  =  PLUS 
RLFA=90 


86 


X — SCRL E  =  5 . OOE  +  Ol  UNITS  INCH. 
Y-5Cfll_E  =  2. OCE+OQ  UNITS  INCH. 


CURRENT  PHASE  :FEM/T=S0LID  :INT.EQ-PLUS 
ZL  =  2 . 0: AL  =  0 . 02 : BL=0. 00 1 : RLFA=90 : NW  =  70 


88 


X-SCAtE=5. 00E-04  UNITS  INCH. 
Y-SCAIE=2.00E*00  UNITS  INCH. 


CURRENT  MRGNITUDE  :FEM/T=SGLID  : I  NT . EQ  =  X 
ZL=2. 0: RL=0. 02: BL=0. 00! : RLFR=60 : NW=70 


89 


X-SCniE.5.00E*01  UNITS  INCH. 

Y-SCBLEsJ.OSE-CO  UNITS  INCH. 

CURRENT  PHRSE  :FEM/T=SOLID  : I  NT . EQ  =  PLUS 
Zl  =  2. 0: RL  =  0. 02: BL  =  0. 001 : RLFR  =  SO:  NW  =  70 


X-SCflLE=5.0CE-0U  UNITS  INCH. 

Y-SCHIE=2.0QE*C0  UNITS  INCH. 

CURRENT  MAGNITUDE  :FEM/T=SOLID  :INT.EQ=X 
ZL=2.0:flL=0.02:8L=0.001:flLFR=30:NW=70 


005  -00U  -003  -002  -001 


X-SCRLE  =  5.  OOE-C'4  UNITS  INCH. 

T-SCRLE= l . OOE+OO  UNITS  INCH. 

CURRENT  MAGNITUDE  :FEM/T=S0LID  :  I  NT . EQ  =  X 
Z  L  = 1 . 0: RL  =  0. 02: B  L  =  0  . 001 : R  L  F  R  =  9  0  : N  W  =  7  0 


93 


X-SCRLE=5 . 0GE+0 1  UNITS  INCH. 

Y-SCRIE=  1 . OOE  +  00  UNITS  INCH. 

CURRENT  PHASE  :FEM/T=SOLID  :INT.EQ=PLUS 
ZL=1.0:AL=0.02:BL=0.001:ALFA=90:NW=70 


94 


X-SCm_E  =  5. 0CE-04 
Y-5CRLE= 1 . OCE+CO 


UNITS  INCH. 
UNITS  INCH. 


CURRENT  MAGNITUDE  :FEM/T=SOLID 
Z  L  =  1 . 0: AL-O. 02: BL=0. 001 : RLFA-SG 


T  V 


95 


unclassified 


naval  POSTGRADUATE  SCHOOL  MONTEREY  CA  f/ 6  20/1R 

rONCEPT  EVALUATION:  Field  FEEDBACK  COMPUTATION  OF  ELECTROMAGNET— EtC tU> 
jUN  AO  B  E  WELCH 


X-SCRIE*5. OOE*01  UNITS  INCH. 

T-5CRtE=l . OOE+OO  UNITS  INCH. 

CURRENT  PHASE  :FEM/T=S0LID  :INT.EQ=PLUS 
ZL=1 . 0: AL=0. 02: BL=0. 001 : ALFA=60: NW=70 


96 


I 


100-  »0-  too-  MO-  too- 


Z 


X-SCfltE*5. 00E-O4  UNITS  INCH. 
Y-5CflLE»1.00E*00  UNITS  INCH. 


CURRENT  MAGNITUDE  :FEM/T=SOLID  :INT.EQ=X 
ZL=1 . 0: AL=0. 02: BL=0. 001 : ALFA=30: NW=70 


97 


X-SCfltE-5.00E*01  UNITS  INCH. 

Y-SCRLE*! .OOE+OO  UNITS  INCH. 

CURRENT  PHASE  :FEM/T=SGLIO  :INT.EQ=PLUS 
ZL=1.0: AL =0.02: BL =0.001 : ALFA= 30: NW= 70 


98 


I 


X-SCBLE-5.006-04  UNITS  INCH. 

T-SCRLE* I . 006*00  UNITS  INCH. 

CURRENT  MRGNITUOE  : FEM/T=SQl I D  :INT.EO=X 
ZL= 1 . 0: AL=0. 02 : BL=0. 001 : RLFR=90 : NW= 100 


X-SCffl.E-S.OOE‘0!  UNITS  INCH. 
Y-SCflLE® 1 . 00E*00  UNITS  INCH. 


CURRENT  PHASE  :FEM/T=SOLID  :INT.EQ=PLUS 
ZL=1.0:AL=0.02:BL=0.001:flLFA=90:NW=100 


.  ■  -  — - + 


100 


X -SCALE »S.OOE -04  UNITS  INCH. 

T-SCALE»I.00E*00  UNITS  INCH. 

CURRENT  MAGNITUDE  : FEM/T=SOL I D  :INT.EQ=X 
ZL=1.0:AL=0.02:8L=0.001:RLFA=60:NH=100 


X-SCPLE»5. OOE+OI  UNITS  INCH. 
Y-SCPLE=I.O0E*O0  UNITS  INCH. 


CURRENT  PHASE  : FEM/T  =  SOl I D  : I  NT . EQ  =  PLUS 
ZL=1.0:flL=0.02:BL=0.001:RLFfi=60:NW=100 


j 


102 


X-SCfllE=5. 00E*01  UNITS  INCH. 
Y-SCfllE=l.OOE*QO  UNITS  INCH. 


CURRENT  PHASE  :FEM/T=SOLID  :INT.EQ=PLUS 
ZL  =  1 . 0: AL  =  0. 02: BL  =  0. 001 : A  L  F  A  =  3  0 : NW=100 


X-SCfil£»5. 00E-04  UNITS  INCH. 
T-SCflLE=5.00E-01  UNITS  INCH. 


CURRENT  MRGNITUOE  :FEM/T=SOLID  : INT.EQ 
ZL=.5:RL=0.01: BL=0. 001: RLFR=90: NW=UO 


111 


mww  wm 


X-SCfiLE»5.00E-*-01  UNITS  INCH. 
Y-SCfllE=5. OOE-O  1  UNITS  INCH. 


CURRENT  PHASE  :FEM/T=SOLID  :INT.EQ=PLUS 
ZL=. 5: AL=0. 01 : BL  =  0. 001  : ALFA=9Q : NW=^0 


106 


-029  -AM  -013  -010  -003  000  005  010  010 


o 

s 


z 


* 
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X-SCflLE=5. 00E-O4  UNITS  INCH. 

Y-SCflLE=5. OOE-Ol  UNITS  INCH. 

CURRENT  MAGNITUDE  :FEM/T=SOLID  :INT.EQ=X 
ZL=. 5: AL=0. 01 : BL=0. 0005: RLF  A  =  90 : NW  =  40 
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-010 


-oos 
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X- SCALE  =  5. OOE *0 1  UNITS  INCH. 
Y-SCAIE=5.00E-01  UNITS  INCH. 


CURRENT  PHASE  :FEM/T=S0LID  :INT.EQ=PLUS 
ZL=. 5: RL=0. 0 1 : B  L  =  0 . 0005: ALFA=90: N  W  =  4  0 


108 


z 


X-SCAIE=5. QOE-04  UNITS  INCH. 

Y-SCAtE=5. 00 E - 0 1  UNITS  INCH. 

CURRENT  MAGNITUDE  :FEM/T=SOLID  :INT.EQ=X 
Zl=.5:AL=0.01:BL=0.0005:RLFR=90:NW=70 


109 


X-SCflLE»5.00E*Ot  UNITS  INCH. 

T-5CRIE=5. OOE-Ol  UNITS  INCH. 

CURRENT  PHASE  :FEM/T=S0LID  :INT.EQ=PLUS 
ZL=. 5:fiL=O.Oi : BL=0. 0005: RLFR=90: NN=70 


-OJS  -OM  -015  -010  -005 


X-SCfiLE=5.00E-OH  UNITS  INCH. 

Y-SCflLE=5.00E-01  UNITS  INCH. 

CURRENT  MAGNITUDE  :FEM/T=SOLID  :INT.EQ=X 
ZL  =  . 5 : AL  =  Q . 0 1 : BL  =  0 . 0005 : ALF A  =  90 : NW= 1 00 


in 


i 


X-SC0LE*5. 00E+01  UNITS  INCH. 

Y-SCfllE=5. 00E-01  UNITS  INCH. 

CURRENT  PHASE  : FEM/T  =  S0L I D .  : I  NT . EQ=PLUS 
ZL=.5:AL=0.01:BL=0.0005:ALFA=90:NW=100 


X -SCALE* 1 . OOE-OH  UNITS  INCH. 

T-5CALES2. OCE -0 1  UNITS  INCH. 

CURRENT  MAGNITUDE  :FEM/T=SOLID  :INT.EQ=X 
Z  L  » . 25: A  l  =  0 . 01 : B  L  =  0  . 00 1 : ALFA  =  90:  N  W  =  4  0 


X-SCBLE-S.OOE‘0!  UNITS  INCH. 


T-SCALE*2.00E-01  UNITS  INCH. 

CURRENT  PHRSE  :  FEM/T  =  SOL I D  :INT.EO=PLUS 
HL=.25:RL*0.01:BL=0.001:RLFR=90:NW=UO 


000-  WO-  MO-  MO-  OIO- 


z 


X-SCfilE-l.OOE-OH  UNITS  INCH. 
T-SCRLE=2.00E-01  UNITS  INCH. 


CURRENT  MAGNITUDE  :FEM/T  =  S0LID  : I  NT . EQ  =  X 
ZL=.25:RL=0.01 :BL =0.0005 :ALFA=9Q :NW= 70 


115 


X-SCfilE«5.00E*01  UNITS  INCH. 
T-SCHLEs2. 00E-01  UNITS  INCH. 


CURRENT  PHR5E  :FEM/T=50LID  :INT.EQ=PLUS 
ZL= . 25 : RL=0 . 0 1 : BL=0 . 0005 : RLFR=90 : NW=70 


116 


X-SCfilE*S.00E+02  UNITS  INCH. 

T-SCRLE=5. 00E-01  UNITS  INCH. 

BOUNDARY  FIELD  MAG.  :  ZL=.5  :  DR=DZ/4 
SOL  I D  =  EXflCT  :  D I AM0NQS=T  :  flL=l/100 


117 


- n 

'  1 1  ’lit  A 


X— SCRLE* 1 . OOE+OO  UNITS  INCH. 

Y-SCRLE=5. 00E-01  UNITS  INCH. 

BOUNDflRT  FIELD  PHR5E  :  ZL=.5  :  DR=DZ/4 
SOL  I D  =  EXRCT  :  TR IRNGLES=T  :  RL=1/100 


118 
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